NASA/CR- 1999-208985 
ICASE Report No. 99-6 



Gas Evolution Dynamics in Godunov-type Schemes 
and Analysis of Numerical Shock Instability 


Kun Xu 

ICASE, Hampton, Virginia 


Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, VA 


Operated by Universities Space Research Association 



National Aeronautics and 
Space Administration 

Prepared for Langley Research Center 
under Contract NAS 1-97046 


Langley Research Center 
Hampton, Virginia 23681-2199 


January 1999 



Available from the following: 


NASA Center for AeroSpace Information (CASI) 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301)621-0390 


National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, VA 22161-2171 
(703) 487-4650 


GAS EVOLUTION DYNAMICS IN GODUNOV-TYPE SCHEMES AND ANALYSIS OF 

NUMERICAL SHOCK INSTABILITY * 


KUN XU+ 

Abstract. In this paper we are going to study the gas evolution dynamics of the exact and approximate 
Riemann solvers, e.g., the Flux Vector Splitting (FVS) and the Flux Difference Splitting (FDS) schemes. 
Since the FVS scheme and the Kinetic Flux Vector Splitting (KFVS) scheme have the same physical mech- 
anism and similar flux function, based on the analysis of the discretized KFVS scheme the weakness and 
advantage of the FVS scheme are clearly observed. The subtle dissipative mechanism of the Godunov method 
in the 2D case is also analyzed, and the physical reason for shock instability, i.e., carbuncle phenomena and 
odd-even decoupling, is presented. 

Key words. Riemann solver, flux vector splitting, flux difference splitting, gas-kinetic scheme, shock 
instability 

Subject classification. Applied Numerical Mathematics 

1. Introduction. In the past decades, tremendous progress has been made in the development of 
numerical methods for compressible flow simulations. Most of them are largely based on the upwinding 
concepts, where the “appropriate” amount of numerical dissipation is implicitly included in the scheme to 
make the smooth shock transition possible [22]. In an earlier paper, we have analyzed the projection dynamics 
in the Godunov method in 1-D case [6, 31], where the implicit viscous term is qualitatively obtained. As a 
continuous effort, in this paper we are going to analyze the dissipative mechanism in the gas evolution stage 
and explain the possible shock instability. 

There arc mainly two kinds of numerical flux functions derived from the inviscid Euler equations. The 
first group is the FVS schemes [24, 26], where the waves generated in the left and right hand side of a cell 
interface move across the cell boundary freely, and the flux has the form F J+ i/ 2 = F Jr (Wj) + 
where Wj and 1 are the left and right states. For example, in the FVS scheme a positive flux function 
from the left state F + (Wj) can move across the cell interface regardless of the existence of the right state 
Wj. j_i. In other words, there is no any dynamical interactions between the left and right states in the FVS 
scheme. Since the flux formulation of the FVS schemes of Steger- Warming and van Leer has the same 
property as that of the Kinetic Flux Vector Splitting (KFVS) scheme based on the collisionless Boltzmann 
equation [8, 14], in this paper the gas evolution mechanism in the KFVS scheme will be analyzed, and 
this analysis can be equally applied to the FVS schemes. People who are not familiar with kinetic schemes 
can find some useful information in [29]. Since the form of the numerical dissipation in the FVS scheme is 
consistent with the Navier-Stokes viscous terms, the robustness of the FVS schemes can be easily understood, 
i.e., the absence of numerical shock instability. However, this advantage also brings the disadvantage for the 
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FVS schemes in solving the physical Navier-Stokes equations, where the large numerical dissipation could 
easily poison the physical viscous term. 

On the other hand, the FDS schemes based on the exact or approximate Riemann solvers, such as the 
Godunov, Roe and Osher methods [6, 16, 21], account for the wave interactions in the gas evolution process. 
Especially for the Godunov method, the exact solution of the Euler equations is used. The wave interaction 
in the FDS scheme can be clearly observed in the Roe’s average W j+i / 2 between the left W 3 and the right 
states Wj+i. In the smooth flow region, there is basically no differences among the Godunov, Roc and Osher 
schemes. For example, the above three schemes can precisely keep a shear layer in the 2D case once the 
shear layer is aligned with the mesh. This fact is consistent with the exact solution of the Euler equations 
Therefore, the FDS schemes can accurately capture the Navier-Stokes solutions in the resolved dissipative 
boundary layer, where the physical viscous term will not be deteriorated from the numerical dissipation. But, 
this advantage also becomes the direct reason for the shock instability, and the so-called odd-even decoupling 
is always associated with the above three schemes [19, 12, 17, 7]. As analyzed in this paper, the artificial 
dissipation of the Godunov method in the 2D case depends not only on the flow distributions, but also on 
the mesh orientation, i.e., in which direction the Riemann solver is applied. Consequently, the form of the 
numerical dissipation in the FDS schemes is basically inconsistent with the Navier-Stokes viscous terms, and 
the Godunov method could generate spurious solutions. 

In this paper, Section 2 is about the physical analysis of the FVS scheme, and Section 3 is related to the 
analysis of the Godunov method. At the same time, the mechanism for the shock instability is presented. 
Section 4 gives reasons for the necessity of developing the hybrid scheme, which combines the advantages of 
both the FVS and FDS methods. The last section is the conclusion and future perspective in the development 
of numerical methods for compressible flow simulations. 

2. Flux Vector Splitting Scheme, In this section, we are going to analyze and understand the 
dynamics in the FVS schemes. The analysis of the KFVS scheme here can be equally applied to understand 
the dissipative mechanism of the Steger- Warming and van Leer methods [24, 26]. 

The flux vector splitting is a technique for achieving upwind bias in numerical flux function, which is 
a natural consequence of regarding a fluid as an ensemble of particles. Since particles can move forward or 
backward, this automatically splits the fluxes of the mass, momentum and energy into forward and backward 
fluxes through the cell interface, i.e., 


F j+ 1/2 “ F + (Wj) + F-(W j+1 ), 


where \Vj represents the mass, momentum and energy densities in the cell j. The equivalence between the 
above splitting mechanism and the collisionless Boltzmann equation was pointed out by Harten, Lax and van 
Leer [8]. Numerically, it is observed that the flux formulation of the van Leer and Steger- Warming schemes is 
almost identical to that of the KFVS scheme [14], In order to clearly understand the dissipative mechanism 
in the FVS schemes, the introduction and analysis of the KFVS scheme is necessary and helpful. A simple 
derivation of the KFVS scheme can be found in [30]. 

It is well-known that the Euler equations can be derived from the Boltzmann equation with the local 
equilibrium distribution function. For an equilibrium state, / is equal to the Maxwellian distribution g , the 
collision term Q(/, /) goes to zero automatically, i.e., Q(g,g) = 0. So, in this case, the Boltzmann equation 
in the 1-D case becomes 

(2-1) ft + uf x = 0, 


2 



where / is the gas distribution function and u is the particle velocity in the ^-direction. Since there is 
no collision term on the right hand side of the above equation, this equation is also called the collisionless 
Boltzmann equation. The KFVS scheme is based on the above governing equation and uses the Maxwellian 
distribution function inside each cell as the initial condition. The KFVS scheme was first derived by Pullin 
and named the Equilibrium Flux Method (EFM) [ 20 ]. With the initial gas distribution function /o(:r, 0 ) at 
t — 0 , the exact solution of the collisionless Boltzmann equation is 

( 2 . 2 ) f = f 0 (x - ut,t). 


With the mass, momentum and energy distribution of a Riemann problem, two equilibrium states at 
x < 0 and x > 0 can be constructed, 


(2.3) 


/ o 


g t , x < 0 

Sr, x > 0 


= Si(l “ H(x)) + g r H(x), 


where H(x) is the Heaviside function. In the 1-D case, the equilibrium state g has the form 
(2.4) g = p (±)*f i e-^ u f+t 2 \ 

7 r 

where the parameters (p, £/, A) are related to the macroscopic mass, momentum and energy densities (p, pU , pe) 
through the relation 

(2-5) p = p, U=U, P - r -—. 

4 pe - ±pU 2 

The parameter K is the dimension of internal variable £ which is a function of the specific heat ratio 7 [29] , 
and £ 2 = £? +£2 + ... +£%. 

With the initial condition in Eq.(2.3), the exact solution of the collisionless Boltzmann equation is 


( 2 . 6 ) f{x,t) = f 0 {x - ut) = gi( 1 - H(x - ut)) + g r H(x - ut). 
Since t > 0 , the above equation can be reformulated as 

(2.7) 


f(x,t) = cfc(l-H(y -u)) + 0 r H(^ - it) 


which is a similarity solution. Based on the solution f(x , t), at a fixed time t we can get the mass, momentum 
and energy distributions in space, 


( 2 . 8 ) 

where d£ = d(,\d^ 2 —d^K and 
(2.9) 


P&’t) . 

pU(x,t) I = I f(x,t)rp a dudfl, 


pe(x, t) 


V’a = 
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For example, for the following initial condition 


(1.0, 0.0, 1.0), ar<0 

(0.1, 0.0, 0.1), z>0 ’ 


( 2 . 10 ) 


(P,U,p) = 


the density distribution at fixed time t is shown in Fig. (5.1), where the collisionless Boltzmann solution 
provides a smooth transition between left and right states. In comparison with the Euler solution, it can 
be clearly seen that the governing equation (2.1) of the KFVS scheme must describe the dissipative flow 
motion instead of the inviscid Euler solution. For the KFVS scheme, the solution f(x,t) at x = 0 is used to 
calculate the fluxes across a cell interface, 


( 2 . 11 ) 



F + (Wj) + F {W j+ i) = f wp„gidud£ + f uil> a g T dud£. 

J u>0 J u<0 


j-f 1/2 


It is true that if the flow remains in a local equilibrium state, Eq.(2.1) is the correct description of the 
inviscid Euler equations. However, even with initial equilibrium states on both sides of a cell interface, the 
collisionless Boltzmann equation cannot keep the equilibrium state at time t ^ 0. For example, the exact 
solution f(x,t) (2.6) is not an equilibrium single Maxwellian at all. Physically, the mechanism for bringing 
the distribution function close to the equilibrium state is the collisions suffered by the molecules of the gas, 
the so-called collision term in the Boltzmann equation. 

In the above analysis, we have only presented the gas evolution model for the KFVS scheme. In order 
to understand the whole numerical solution, we have to include the dynamics from the projection stage too. 
Physically, the projection stage is mainly to use the total mass, momentum and energy inside each cell to 
reconstruct an equilibrium state, which is equivalently to introduce particle collisions in the gas system to 
make the transition from the non-equilibrium state f(x,t) to an equilibrium state g possible. Therefore, the 
projection mechanism can be precisely described by the following equation 


( 2 * 12 ) f t = Q {f,f), 

where stiff collision term Q(/, /) make the transition from f to g instantaneously. 

Combining the gas evolution model (2.1) and the projection model (2.12), the real governing equation 
of the KFVS scheme can be physically described as the modified “BGK” equation [1], 


(2.13) 


df 

dt 


+ u 


df 

dx 


9~ f 
At ’ 


where the real particle collision time r in the BGK model is replaced by the numerical time step At. In should 
be emphasized again that the above equation can be equally applied to understand the Steger- Warming and 
van Leer’s FVS schemes. After Eq.(2.13) is obtained, all theoretical results related to the BGK model of the 
approximate Boltzmann equation can be applied here [1]. For example, the modified BGK model have the 
intrinsic viscosity coefficient rj and heat conductivity coefficient /c, 


(2.14) 


g = pAt, 


and 

(2.15) K =E±AJL p At, 

2 rn 
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where p is the pressure, k the Boltzmann constant, and m the molecule mass. More importantly, even 
though the magnitude of the viscosity and heat-conductivity coefficients are always proportional to the time 
step, the format of the numerical viscous terms from (2.13) are always consistent with those in the Navier- 
Stokes equations. The artificial dissipation exists in both smooth and discontinuity flow regions, and it is 
independent of the mesh-orientation. As a result, the FVS schemes totally avoid the carbuncle phenomena 
and odd-even decoupling [7]. 

As analyzed in [29], the fact p ~ At cannot be changed by increasing the order of the scheme, it is from 
the free transport dynamics. Although the FVS schemes arc robust, the accuracy of the schemes is limited. 
The large artificial dissipation in the FVS scheme could easily deteriorate the physical viscous term in the 
Navier-Stokes flow calculations [13, 27, 11], 

3. Flux Difference Splitting Scheme. The well-known FDS schemes include the Godunov, Roc and 
Osher’s methods. The analysis of the Godunov method in this section can also be applied to the Roc and 
Oshcr schemes in the explanation of the shock instability. 


3.1. The Godunov Gas Evolution Model. The Riemann problem is defined as an initial value 
problem for the Euler equations. In the 1-D case, with the following initial condition at t = 0, 


(3.1) 


Mp) 


(PL,U L ,p L ), x < 0, 
(Pr,Ur,pr), x > 0, 


the entropy-satisfying solutions arc the following: the left state {pl,Ul,Pl) is connected to the right state 
(Pr,Ur,Pr) by a 1-shock or 1-rarefaction wave, a 2-contact discontinuity, and a 3-shock or a 3-rarcfaction 
wave. The 2-contact discontinuity separates two constant states (p/, U*,p*) and (p//, £/*, p*), so that (17, p) 
are continuous across the contact discontinuity. For example, in Fig.(5.2), the 1-wavc is a rarefaction and 
the 3- wave a shock. There is standard technique to obtain the solutions around a contact discontinuity [25]. 
Basically, the Godunov method is closely following the solution of the inviscid Euler equations. 

Once the solution of the Riemann problem is obtained, the flow variables at x — 0, i.c., ( pj,U*,p *) in 
Fig. (5.2), are used to construct fluxes. The Godunov method uses these fluxes across each cell interface to 
update the flow variables. In order to understand the gas evolution model in the Godunov method, we have 
to understand the underlying assumption in the process of evaluating the fluxes. For the Euler equations, 
the corresponding fluxes are obviously 



f F P \ 

f PU \ 

(3.2) 

FpU 1 = 

pU 2 4 p 


\ F pJ 

U PU 3 + ^pu) 


For example, in Fig. (5. 2), the fluxes at x = 0 can be written as (pjU*, piU* 2 + p*, |p/7* 3 + ^yp*t/*) T . 
From the flow variables (p,U,p) to the above flux functions (F p , F pU , F pe ), it is implicitly assumed that 
the gas is staying in an equilibrium state. In other words, whatever the real flow situations is, the state 
(pi,U*,p*) at the cell interface always corresponds to a single Maxwellian distribution function g in the 
Godunov method, even for the fluid inside a numerical shock layer, see Fig. (5.3). On the contrary, the 
KFVS scheme uses two half Maxwcllians gi\ u >o and g r \ u < o in the flux evaluation (2.11). 

As we know, the intermediate points in the numerical shock region have to be regarded as points inside 
a shock structure, where the dissipative mechanism is crucially important to translate kinetic energy into 
thermal energy and make the smooth shock transition. But, the gas evolution model of the Godunov method 
inappropriately assigns the equilibrium state there at the cell interfaces. In the special case, if the shock 
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is located exactly at a cell interface, the Godunov method can precisely capture a thickless shock, such as 
that shown in Fig. (5. 4). However, if we look at the Riemann solution in this case, we cannot even find 
an unique state (p,t/,p) at the cell interface. So, in this case the Godunov method uses the continuous 
flux function directly instead of constructing any equilibrium state, subsequently avoids mistakes. In the 
general case, once the numerical shock covers a few grid points or the shock is unsteady, the weakness of the 
Godunov method automatically appears. The ignorance of the dissipation in the gas evolution stage in the 
discontinuous region is one of the disadvantages of the Godunov method. Fortunately, as analyzed in [31], 
in the 1-D case the projection stage provides the dissipation needed in the shock region for the Godunov 
method to capture numerical shocks. 

3.2. The Godunov Projection Dissipation In 2D Case. In the 2-D case, for any shock wave 
with upstream and downstream velocity (Ui(x, y), U 2 (a:, y)) and density distributions (p^rr, y), p 2 (x, y)), 
the dissipation provided in the projection stage inside any finite elemet AV is 

(3.3) A E k =E k - E k) 

where E k is the total kinetic energy before the averaging, E k is the kinetic energy after the averaging, and 
the difference is the kinetic energy lost during the averaging process in that volume, which is translated into 
the thermal ones. Suppose looked in the n-dircction there is a flow velocity difference inside a small element, 
such as that shown in Fig. (5. 5). For simplicity, the volumes occupied by (pi,Ui) and (p 2 ,U 2 ) are assumed 
to be the same and equal to 1/2. From the conservations of the mass and momentum, we have 

P = “(pi + P 2 ) 

and 

pU = 2 (Pi^i +P2U2). 

As a result, the kinetic energy before the averaging is 

E k — ~(piUi 2 + P2U2 2 ), 

and after the averaging 

E k = IpU 2 . 

With the above relations, the difference in the kinetic energy can be obtained, 

4 Pi + P 2 

where n is the direction in which the fluid velocity distributions are appeared to be Ui , U 2 . So, the projection 
stage in the Godunov method provides the follwoing subtle dissipation 

(3.4) dissipation ~ P 1 ? 2 m 1 — U 2 )l. 

P1P2 

Again, h is also the direction in which the Riemann solution is evaluated. For the 1-D flow, (U x - U 2 ) is on 
the same direction as n and the dissipation is always being proportional to (Ui — U 2 ) 2 • So, in the 1-D case, 
the velocity difference provides intrinsic dissipation to the Godunov method in the shock region. However, in 
the 2-D case, the velocity distributions Ui and U 2 depend on which direction you look at the flow, and the 
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dissipation in the Godunov method is mesh-oriented. As shown in Fig. (5.6), for the same fluid distribution 
in the 2-D case, i.e., a normal shock in the ^-direction, the velocity distribution (Ui - U 2 )^ is different 
when looked in the different directions n. So, the dissipative term in Eq.(3.4) is directionally dependent* 
Therefore, the form of the dissipation provided in the projection stage in the Godunov method cannot be 
consistent with the Navier-Stokes dissipation at all. As shown in the case (b) of Fig. (5.6), for a 2-D normal 
shock in the ^-direction, the dissipation provided in the y-direction becomes (Ui - U 2 )y ~ 0, because U x 
and U 2 appear as the parallel flows in that direction. 

3.3. The Explanation of Shock Instability. The documented observations of the shock instability 
arc scattered in many literatures. The influential paper by Quirk for the first time systematically presented 
the observation and analysis [19] . After that, the explanation of this numerical shock behavior has attracted 
much attentions. A short list of references includes [4, 7, 12, 17, 29, 23], and references therein. 

In order to have the shock instability, such as the carbuncle phenomena or odd-even decoupling, the 
numerical fluid has the following properties [19, 12]; 

(1) . Shock propagates and is aligned with the mesh. 

(2) . As pointed out by Liou [12], the mass flux due to a pressure difference is not zero, such as D ^ ^ 0 in 

the dissipative term V = A p + U + D^Ap of the full flux function F J+l/ 2 = \{F j + F j+1 ) + V. 

(3) . It is observed that the shock instability first happens in the downstream region, sec Fig. (5) in Quirk’s 
paper [19], then amplifies rapidly to the whole shock region. The formation of the instability is independent 
of the time step used in the calculation. 

From the point (1), if the shock front is aligned with the mesh, the projection dissipation (3.4) can be 
provided in the direction normal to the shock front only, such as in the ^-direction in Fig. (5.6). There is 
almost no projection dissipation provided in the direction along the shock front due to the equal velocity. 

Now suppose that a stationary shock is staying in the ^-direction. Inside the shock region, due to 
numerical perturbations, it is possible that certain fluid element will move from cross section Sa to S& 
through the streamline converging or diverging, such as that shown in Fig.(5.7). For this element, the 
momentum conservation gives 

(3.5) pS A + pU 2 S A + pdS A = (p + dp)(S A + dS A ) + (p + dp)(U + dU) 2 (S A + dS A ), 

where Sb = Sa + cISa , and the continuity equation d(pUSA ) = 0 has 

pU (ISa + pSacIU + SaUcLp = 0. 

Multiplying the above equation by U and Subtracting it from Eq.(3.5), we have 

dp — — pUdU , 

which means that the pressure change is always out of phase with the velocity change. From aerodynamic 
theory for the qusi-one- dimensional flow, we then have 

Case(l). For 0 < M < 1 (downstream side), an increase in velocity (positive dU) is associated with a 
decrease in area (Sa), and vice versa. Therefore, the velocity increases (pressure decreases) in a converging 
streamline case and decreases (pressure increases) in a diverging streamline case. 

Case(2). For M > 1 (upstream side), the velocity increases (pressure decreases) in diverging streamline case 
and decreases (pressure increases) in the converging streamline case. 

Once this kind of perturbation exist inside the numerical shock front, in the subsonic side (close to 
downstream), an increase of velocity is associated to a decrease of pressure in the converging streamline case, 
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see Fig. (5.8). Physically, the shear viscosity will strongly take effect to reduce the velocity differences and 
stablize the numerical shock structure in this situation. However, numerically, following the inviscid Euler 
solution in the y-dircction (no projection dissipation provided in this direction) will give a different answer. 
The exact 1-D Ricmann solver in the y-direction does not recognize the x-component velocity differences 
due to the inviscid nature. The only dynamical influence involved in the Ricmann solver in this direction is 
the pressure difference. Since the fluid in the surrounding cells has a higher pressure, it will move toward 
the center. As a consequence, the fluid in the central cells passes through an even narrower cross section, 
and the fluid speed gets even higher due to the converging of streamlines. At the same time, the pressure 
in the central cells is even lower. So, the use of the exact Riemann solver in the y-direction could only help 
to amplify the initial perturbations, and the numerical fluid will speed up and penetrate the shock layer in 
the subsonic side first to form instability. Numerically, it does observed that the instability happens first in 
the downstream side [19]. With the converging of streamlines, the size of the cross section becomes small 
and will reach the cell size eventually to form the saw-tooth like patterns. In the supersonic side of the 
numerical shock structure, such as the left part of Fig. (5. 7), in case of the converging of streamlines, an 
increase in velocity will be associated with an increase in pressure. So, the central high pressure gas will 
expand and stop the converging of the streamlines. Therefore, in the upstream side M > 1, the structure 
is theoretically stable to small perturbations. However, once the perturbation in the subsonic region gets 
amplified (saw-tooth like profile), the flow in the supersonic side can hardly match the downstream flow 
variations. Then, the shock strength becomes different along the y-dircction. As a result, the different 
shock propagating speed will eventually destroy the whole numerical shock structure. For the odd-even 
decoupling case, the streamline can diverge or converge, and the instability happens in both cases. For the 
carbuncle phenomena, the streamline after the normal shock in front of the cylinder basically diverges, So, 
the formation of instability in this case can be more difficult than that in the odd-even decoupling case. But, 
physically the mechanism of the shock instability is basically the same. The above analysis also validates the 
numerical observation that the shock instability does not happen in the flow calculation with the unstructurc 
mesh, where the shock front can hardly be systematically aligned with the irregular mesh 1 . In this case, the 
projection dissipation exist all the time. Also, if the mesh is not aligned with the shock front, as shown in 
the cases of (c) and (d) in Fig. (5. 6), the projection dissipation is provided in both directions to eliminate 
the shock instability. 

The above explanation validates Liou’s conjecture that T> ^ ^ 0 is accompanied by an instability [12]. 
If = 0, the initial pressure difference in the y- direction will not push the fluid toward the central cells, 
and no instability will be formed. Unfortunately, even for the exact Riemann solver, T> ^ ^ 0 always holds. 
So, the shock instability is intrinsically rooted in the Godunov method if the inviscid Euler equations are 
considered as the governing equations in the numerical shock region in the y-direction. The above explanation 
has no conflict with the stability of inviscid shock [10], where a discontinuity shock (without structure) is 
perturbed and is stable. Here, in our analysis we are concerning with the mis-use of the governing equations 
in the dissipative shock region and the reasons for the shock structure instability. For the FVS scheme, the 
story is different. Once there is any perturbation in the velocity difference, the strong shear viscous term 
will take effect and suppress the amplification of the instability. 

For the shock instability case, Quirk used a saw-tooth initial condition to test the response from the 
Roe’s Riemann solver [19]. It was observed that the initial density perturbation goes to p°° = /5°/a 2 , where 

1 Philosophically, the randomness of the unstructure mesh is more consistent with the homogeneous physical space (no 

preferred direction) than the well-organized rectangular mesh, where the homogeneous physical space is systematically distorted. 
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a is the local sound speed. At the same time, the pressure perturbation gets decaying. Similar analysis 
was also given in [17]. In these analysis, in order to get the shock instability, as pointed out by Quirk [19], 
a systematic pressure perturbation is needed to amplify p°°. Then, the question is where the systematic 
perturbation comes from. Actually, Quirk’s analysis is only based on the solution of an initially saw-tooth 
shear layer, it totally ignores the basic fact that the instability happens in the numerical shock region, where 
the ^-direction fluid velocities change dramatically across the shock layer. Therefore, the shear instability 
analysis can hardly explain the shock instability. 

In conclusion, the odd-even decoupling and the carbuncle phenomena arc due to the numerical dis- 
cretization of the inviscid Euler equations in the shock region and the special dissipative mechanism of the 
Godunov method in the projection stage. Condition (1) is necessary to have this kind of phenomena to take 
place. In order to stabilize the shock front in y-direction, a direct way is to add dissipation through entropy 
fix, such as to solve the viscous shear layer equation. However, as realized by Quirk[19], the entropy fix here 
is added to the linear wave rather than the nonlinear ones. So, it is more precisely to say that we need to 
solve viscous governing equations once the “appropriate” dissipation from the upwinding is not enough. 

4. The Necessity of Hybrid Scheme. As analyzed in the previous sections, the dissipative mecha- 
nism in the FVS and FDS schemes are different. Due to the free transport in the gas evolution stage and the 
projection dissipation, the FVS scheme is basically solving the Navier-Stokes-type equations and the dissipa- 
tion is present in both the smooth and discontinuity flow regions. Consequently, the viscous terms help the 
FVS scheme avoid any shock instability. On the contrary, for the Godunov method, since the exact inviscid 
Euler solution is used in the gas evolution stage, the dissipation is mainly contributed from the projection 
stage. However, the projection dissipation is not only a function of flow variations, but also a function of 
mesh orientation. This fact makes the numerical dissipation of the Godunov method be inconsistent with the 
Navier-Stokes terms. In other words, the dissipation strongly depend on the direction in which the Riemunn 
solver is applied. For example, even inside the numerical shock region, if the mesh is aligned with the shock 
front, in the direction along the shock front the Godunov method interprets the dissipative shock layer as 
an inviscid shear layer. As a result, the shock instability emerges. 

In order to construct a scheme with both robustness and accuracy, a recent trend in the development of 
upwind schemes has centered around the construction of hybrid flux-splitting formulations which combine 
the accuracy of the FDS approach in the resolution of shear layer and the robustness of the FVS method in 
the capturing of strong discontinuity [3, 4, 11, 13, 15, 2, 9]. For example, starting from the FVS scheme, in 
order to reduce the free penetration mechanism, the correlation and collision between left and right moving 
waves have to be introduced. The constructions of the values M 1/2 ,^i /2 at a cell interface in AUSM-type 
schemes arc largely based on this physical consideration^ 1, 12]. Moschetta and Pullin [15] tried to use 
Osher s numerical flux to replace linearly degenerate wave in the KFVS scheme. Unfortunately, this fix 
introduces the shock instability [7]. This fact validates our analysis in the last section again. From the FDS 
scheme, in order to eliminate the shock instability, a direct remedy is to introduce additional viscous term. 
As a routine practice, additional dissipation had been introduced to stabilize discontinuous solution in the 
early development of high resolution scheme [28]. The dynamical effect of using two waves instead of three 
waves in the HLL scheme is basically to smear some waves and introduces dissipations in the gas evolution 
stage [8]. The same philosophy is used in Marquina’s flux function [3], where Steger and Warming’s FVS flux 
is modified. A good hybrid scheme depends on a smart weight function to identify where the flow is smooth 
or where the flow is discontinuous. The weakness for these hybrid methods is that they are not based on 
any rigorous governing equations, rather than ad hoc fixes whenever unsatisfactory solution appears. Also, 
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no uniform criteria have ever been established so far for the optimum choice [19], Numerically, as tested 
by Grcssier and Moschetta [7], almost no hybrid schemes [12, 2, 15] could really remove the pathological 
behavior of upwind schemes in the refining mesh study if the scheme does not use the FVS-typc homogeneous 
dissipative term. 

Besides the above approaches, in the gas evolution stage the gas kinetic BGK scheme is based on the 
collisional Boltzmann equation [29, 32], 


df 

dt 


+ u 


df 

dx 


9~ f 

? 


where the viscosity and heat conductivity coefficients in the BGK scheme are related to the particle collision 
time r, instead of the time step At in the KFVS method. So, in the smooth region, accurate Navier-Stokes 
solutions can be obtained by the BGK method by appropriately choosing the collision time r according to the 
physical viscosity coefficient [29]. For the Euler solution, the small but non-zero r provides a homogeneous 
background dissipation, which subsequently avoids the shock instability. The resolution of the BGK solution 
is as good as those from the best FDS method. In the discontinuity region, the BGK method uses large 
artificial particle collision time r ^ At, the numerical dissipation is obtained to construct a stable shock 
structure, such as the KFVS scheme. Therefore, the BGK method naturally provides a genuine nonlinear 
combination of the FVS or FDS schemes [29]. Different from any attempt in the modification and mixing 
of the FVS and FDS flux functions, the BGK method is based on the a real physical governing equation. 
Whatever the value of the particle collision time is used, the format of the viscous terms is always consistent 
with the Navier-Stokes equations. It is basically the direct reason for the robustness and accuracy of the 
BGK method. Even though there is a “tunable” parameter t here, physically the fluid does need us to use 
different governing equations in the smooth and discontinuity regions. 


5. Conclusion. Quite often, the requirements for the robustness and accuracy in the design of a 
numerical scheme are in conflict with each other: if a scheme is robust, it is unnecessarily diffusive; if a 
scheme is accurate, it loses robustness. It is well-known that the inviscid Euler equations cannot give a 
correct representation of fluid motion in the discontinuity flow region. In order to capture the discontinuity 
solution appropriately, the explicit dissipation has to be added in the gas evolution stage. For the FVS 
schemes, due to the free transport gas evolution model, the numerical viscous terms of the scheme arc 
consistent with the Navier-Stokes equations and the dissipation appears in both smooth and discontinuity 
regions. For the Godunov method, the numerical dissipation in the multidimensional case is mesh-oriented, 
and it is not consistent with the Navier-Stokes equations. As a result, in the 2-D case, even inside the shock 
layer, the numerical dissipation is reduced to the minimum level when the shock front is aligned with the 
mesh, and the shock instability forms. 

In order to develop more accurate and robust numerical schemes for the compressible flow simulation, 
first we have to have a physically consistent governing equation. Specifically, we have to solve the viscous 
governing equations directly in the discontinuity region. No perfect scheme, or the Riemann solver, will be 
obtained if only the Euler equations are regarded as the governing equations. It is naive to believe that 
numerics only could provide the “appropriate” and “consistent” physics once it is necessary. For the gas- 
kinetic BGK scheme, we always solve the viscous governing equations and the dissipation is controlled by the 
collision time. With the increasing of computer power and the continuing refinement of the numerical mesh 
size, the numerical error will become less and less, and the consistency of the digital physics in the scheme 
with the real physics will become more and more important, especially for the complicated flow calculations. 
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Equilibrium states at cell interfaces 



Dissipative Numerical Shock Region 


Fig. 5.3. Godunov Gas Evolution Model. 


Non-uniquesness states at a cell interface 



Fig. 5.4. A stationary shock front located exactly at a cell interface. 
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Fig. 5.7. A normal shock in the redirection in the 2-D case 
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Fig. 5.8. Schematic explanation of shock instability 
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